Identifying the impacts of extreme weather

A case study of Houston, Texas

Link to GitHub repository: https://github.com/aakriti-poudel-chhetri/impacts-of-extreme-weather

Introduction

Climate change is increasing the frequency and intensity of extreme weather events, with devastating impacts. In February 2021, Texas experienced a major power crisis caused by three severe winter storms sweeping across the United States on February 10–11, 13–17, and 15–20. These storms led to widespread electricity outages, particularly in the Houston metropolitan area, but the spatial and socioeconomic distribution of these outages remains uncertain.

This analysis seeks to understand the spatial and socioeconomic impacts of the blackout by estimating the number of homes that lost power and examining whether these impacts were disproportionately felt across communities. Using remotely-sensed night-lights data from the VIIRS VNP46A1 satellite, areas of lost electric power are detected by comparing night-light intensity before and after the storm. These blackout areas are then linked with OpenStreetMap building data to estimate affected homes and with US Census Bureau data to explore potential socioeconomic patterns, such as differences in median household income between affected and unaffected census tracts.

The study focuses on several related questions such as how many homes likely lost power during the winter storms; which geographic areas experienced measurable reductions in night-light intensity once highways are excluded; which census tracts were affected; and how does blackout exposure relate to median household income. Together, these questions aim to determine whether the blackout disproportionately impacted higher or lower income neighborhoods, and identify where, and to what extent, communities across Houston experienced power outages.

This approach highlights both the spatial distribution of storm impacts and the limitations of using remote-sensing data for assessing power outages.

Set up

Let’s start by importing all the necessary libraries.

library(tidyverse)
library(sf)
library(tmap)
library(terra)
library(stars)

About data

The data for this study were retrieved from multiple sources listed below.

Night lights - The NASA’s Worldview is explored for extracting the data around the day of the storm. There are several days with too much cloud cover to be useful, but 2021-02-07 and 2021-02-16 provides two clear, contrasting images to visualize the extent of the power outage in Texas.

VIIRS data is distributed through NASA’s Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC). Many NASA Earth data products are distributed in 10x10 degree tiles in sinusoidal equal-area projection. Tiles are identified by their horizontal and vertical position in the grid. Houston lies on the border of tiles h08v05 and h08v06. These two tiles per date were pre-downloaded and provided by the team.

Roads - We used Geofabrik’s download sites to retrieve a shapefile of all highways in Texas and prepared a Geopackage (.gpkg file) containing just the subset of roads that intersect the Houston metropolitan area.

Houses - The data were downloaded from Geofabrik and processed to create a GeoPackage containing only residential buildings within the Houston metropolitan area.

Socioeconomic - The socioecoomic data were obtained from the U.S. Census Bureau’s American Community Survey for census tracts in 2019.

Let’s read the datsets.

# Read night lights data
tile05_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v05.001.2021039064328.tif'))
tile06_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v06.001.2021039064329.tif'))
tile05_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather','data', 'VNP46A1', 'VNP46A1.A2021047.h08v05.001.2021048091106.tif'))
tile06_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021047.h08v06.001.2021048091105.tif'))


# Read roads data
highways <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'gis_osm_roads_free_1.gpkg'),
                 query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass = 'motorway'") %>%
  st_transform(crs = 'EPSG:3083')


# Read houses data
houses <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data',
                          'gis_osm_buildings_a_free_1.gpkg'),
  query = "SELECT *
    FROM gis_osm_buildings_a_free_1
    WHERE (type IS NULL AND name IS NULL)
      OR type IN ('residential', 'apartments', 'house', 'static_caravan', 'detached')") %>% 
    st_transform('EPSG:3083')


# Explore the contents of the geodatabase
socioeconomic <- st_layers(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'))

# Extract the census tract information
census_tract <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'ACS_2019_5YR_TRACT_48_TEXAS') %>% 
  st_transform('EPSG:3083')

# Extract the income information
income <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'X19_INCOME')

First, the raster tiles were merged and cropped to the Houston bounding box to focus the analysis on the study area and ensure that only valid night-light data within Houston were used for assessing power outages.

# Merge two raster data of February 7
nightlight_07 <- st_mosaic(tile05_07, tile06_07)

# Merge two raster data of February 16
nightlight_16 <- st_mosaic(tile05_16, tile06_16)

# Create Houston bounding box
houston_bbox <- st_bbox(c(xmin = -96.5, xmax = -94.5, ymin = 29, ymax = 30.5),
                        crs = 'EPSG:4326')

# Crop all raster data to Houston area
nightlight_07_crop <- st_crop(nightlight_07, houston_bbox)
nightlight_07_crop[(nightlight_07_crop > 1000) | (nightlight_07_crop <= 0)] <- NA

nightlight_16_crop <- st_crop(nightlight_16, houston_bbox)
nightlight_16_crop[(nightlight_16_crop > 1000) | (nightlight_16_crop <= 0)] <- NA

Here, we visualize the comparison between night-light intensities in Houston before and after the February 2021 winter storm. By mapping the February 7 (pre-storm) and February 16 (post-storm) night-light data side by side, the analysis shows areas where lighting and likely power availability changed, helping to identify regions affected by the blackout.

# Map for February 7 (before storm)
beforestorm <- tm_shape(nightlight_07_crop) +
  tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
            col.legend = tm_legend(title = 'Night light\nintensity',
                                   position = tm_pos_out('right'))) +
  tm_title('Houston Night Lights\nFebruary 07, 2021 (Before storm)',
           size = 1.3) +
  tm_compass(position = c('left', 'top'),
             type = 'arrow',
             size = 3,
             text.size = 0.6,
             color.dark = 'grey50',
             text.color = 'white') +  
  tm_scalebar(position = c('right', 'bottom'),
              text.size = 0.8,
              color.dark = 'grey50',
              text.color = 'white') +  
  tm_layout(bg.color = 'white',           
            outer.bg.color = 'white',      
            frame = FALSE,
            legend.show = FALSE)


# Map for February 16 (after storm)
afterstorm <- tm_shape(nightlight_16_crop) +
  tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
            col.legend = tm_legend(title = 'Night light\nintensity',
                                   position = tm_pos_out('right'))) +
  tm_title('Houston Night Lights\nFebruary 16, 2021 (After storm)',
           size = 1.3) +
  tm_compass(position = c('left', 'top'),
             type = 'arrow',
             size = 3,
             text.size = 0.6,
             color.dark = 'grey50',
             text.color = 'white') +  
  tm_scalebar(position = c('right', 'bottom'),
              text.size = 0.8,
              color.dark = 'grey50',
              text.color = 'white') +  
  tm_layout(frame = FALSE,
            legend.show = FALSE)

legend_map <- tm_shape(nightlight_16_crop) +
  tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
            col.legend = tm_legend(title = 'Night light\nintensity',
                                   position = tm_pos_out('right'))) +
  tm_layout(legend.only = TRUE, legend.outside = TRUE, legend.position = c("right", "center"))

tmap_arrange(beforestorm, afterstorm, legend_map, ncol = 3)

We will now create a blackout mask by identifying areas of Houston where night light intensity decreased significantly after the storm, indicating potential power outages. By calculating the difference between before and after night light, cropping to the study area, filtering out minor changes, and vectorizing the result, the mask provides a spatial representation of blackout affected regions for further analysis.

# Calculate the difference in light intensity
nightlight_diff <- nightlight_07 - nightlight_16

# Crop and reclassify the raster data 
mask <- st_crop(nightlight_diff, houston_bbox)
mask[mask < 200] <- NA

# Vectorize the blackout mask
blackout <- st_as_sf(mask,
                     as_points = FALSE,
                     merge = TRUE) %>%
  st_make_valid() %>%
  st_transform(crs = 'EPSG:3083')

In this step, highway geometries are combined and buffered to exclude areas where changes in night light intensity could be due to traffic rather than power outages. By removing blackout areas within 200 meters of highways, the analysis isolates regions where lighting loss more likely reflects true residential or commercial power outages.

# Combine all highway geometries into one
highways_union <- st_union(highways)

# Create 200m buffer around all highways
highways_buffer <- st_buffer(highways_union, dist = 200)

# Find areas that experienced blackouts further than 200m from highways
blackout_far <- st_difference(blackout, highways_buffer)

# Visualize the blackout in Houston
tm_basemap('OpenStreetMap') +
  tm_shape(blackout_far) +
    tm_polygons(fill = 'ivory', fill_alpha = 0.8, col = 'black') +
    tm_title('Blackout areas in Houston', fontface = 'bold', size = 1) +
    tm_layout(legend.show = TRUE) +
    tm_scalebar(position = c('left', 'bottom')) +
    tm_compass(position = c('right', 'top'))

Now estimate the number of homes in Houston that lost power during the February 2021 winter storms. By identifying which building footprints intersect with the blackout mask, the analysis isolates the homes likely affected by outages. Converting buildings to centroids allows for precise spatial calculations and comparisons between impacted and non-impacted houses, resulting in an estimated total of approximately 157,970 homes without power.

# Keep buildings that intersect with blackout areas
houses_blackout <- st_intersects(houses, blackout_far)

# Check if each building is in a blackout area
in_blackout <- lengths(houses_blackout) > 0

# Filter buildings in blackout areas
blackout_houses <- houses[in_blackout, ]

# Number of impacted buildings
n_blackout_houses <- nrow(blackout_houses)

# Convert all buildings to centroids
houses_points <- houses %>% 
  st_centroid()

Warning: st_centroid assumes attributes are constant over geometries

# Houses impacted by blackouts
houses_impacted <- blackout_houses %>% 
  st_centroid()

Warning: st_centroid assumes attributes are constant over geometries

# Houses not impacted by blackouts
houses_not_impacted <- houses_points %>% 
  filter(!osm_id %in% blackout_houses$osm_id)

This map visualizes the spatial distribution of homes in Houston that lost power during the February 2021 winter storms. Impacted houses are highlighted in red, while unaffected houses are shown in blue, allowing for a clear comparison of affected versus unaffected areas across the city. The map provides a detailed spatial perspective on the extent of the blackout and helps identify which neighborhoods were most impacted.

# Create the map
tm_basemap('OpenStreetMap') +
  tm_shape(houses_not_impacted) +
    tm_dots(fill = 'midnightblue', size = 0.02, fill_alpha = 0.8) +
  tm_shape(houses_impacted) +
    tm_dots(fill = 'firebrick', size = 0.02, fill_alpha = 0.8) + 
  tm_title('Houses in Houston that lost power',
           position = tm_pos_out('center', 'top'), 
           fontface = 'bold', size = 1.2) +
  tm_layout(legend.outside = TRUE,
            legend.outside.position = 'right',
            legend.title.size = 1,
            legend.text.size = 0.8) +
  tm_add_legend(labels = c('Impacted houses', 'Not impacted houses'),
                fill = c('firebrick', 'midnightblue'),
                col = c('firebrick', 'midnightblue'),
                alpha = c(0.2, 0.7, 0.5),
                title = 'Legend') +
  tm_scalebar(position = c('left', 'bottom'),
              text.size = 0.6) +
  tm_compass(position = c('right', 'bottom'),
             type = '4star',
             size = 2)

The following process links census tract boundaries with median household income data to analyze the socioeconomic characteristics of neighborhoods affected by the blackout. By joining the income dataset to the census tract geometries, the analysis enables comparisons between impacted and non-impacted areas and supports investigation of whether power outages disproportionately affected higher- or lower-income communities.

# Rename the column in income dataset to match with the census tract
income_renamed <- income %>% 
  rename("GEOID_Data" = "GEOID")

# Left join the two datasets by the geometry
census_income <- census_tract %>% left_join(income_renamed, by = 'GEOID_Data')

Let’s verify CRS. Verifying CRS alignment is essential for accurate spatial analysis, as mismatched projections can lead to incorrect geometry intersections, area calculations, or mapping errors.

# Check CRS of the dataset
if(st_crs(census_income) == st_crs(houston_bbox)){
print("Coordinate reference systems match!")
} else {
warning("Update coordinate reference systems to match!")
}

Warning: Update coordinate reference systems to match!

By aligning the coordinate systems and cropping to the Houston study area, census tracts and blackout affected houses are prepared for spatial comparison. Intersections between the two datasets identify which tracts experienced power outages, and a new column flags these impacted areas. Filtering to only the affected tracts enables a targeted analysis and visualization of how the blackout influenced neighborhoods across the city.

# Transform the CRS and crop to houston bounding box
census_income <- st_transform(census_income, crs = 'epsg:4326') 

# Crop census income to Houston
census_income_crop <- census_income %>% 
  st_crop(houston_bbox)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

# Transform and clean blackout houses
blackout_houses_transformed <- blackout_houses %>% 
  st_transform(crs = 'epsg:4326')

# Identify the blackout areas
census_blackout <- st_intersects(census_income_crop, blackout_houses_transformed)

# Create the column
census_blackout_col <- census_income_crop %>% 
  mutate(blackout = lengths(census_blackout) > 0)

# Select only the values that are true for plotting
census_blackout_col_filtered <- census_blackout_col[census_blackout_col$blackout, ]

The following map highlights Houston census tracts that experienced power outages during the February 2021 winter storms. Affected tracts are shown in red, providing a clear visualization of the spatial distribution of blackouts across the city. This allows for easy identification of neighborhoods impacted by the storm and supports further analysis of socioeconomic patterns in relation to power loss.

# Map the census tracts that lost power
tm_basemap('CartoDB.VoyagerNoLabels') +
  tm_shape(census_blackout_col_filtered) +
  tm_polygons(fill = 'blackout',
              fill.scale = tm_scale_categorical(
                values = c('TRUE' = 'firebrick'),
                labels = c('TRUE' = 'Blackout')),
              fill.legend = tm_legend(title = 'Blackout status'),
              col = 'grey50',
              lwd = 0.2) +
  tm_title('Census tracts in Houston that lost power',
           position = tm_pos_out('center', 'top'),
           fontface = 'bold',
           size = 1) +
  tm_layout(legend.outside = TRUE,
            legend.outside.position = 'right',
            legend.bg.color = 'white',
            legend.bg.alpha = 0.8) +
  tm_compass(position = c('right', 'bottom'),
             type = 'arrow',
             size = 1) +
  tm_scalebar(position = c('left', 'top'),
              text.size = 0.8)

This plot compares the distributions of median household income between Houston census tracts that experienced blackouts and those that did not during the February 2021 winter storms. By visualizing income differences, it allows for assessment of whether power outages disproportionately affected higher- or lower-income neighborhoods and provides insight into potential socioeconomic disparities in storm impacts.

ggplot(census_blackout_col %>% filter(!is.na(B19013e1)), # Remove NA values
                                      aes( x = blackout, y = B19013e1, fill = blackout)) +
  geom_boxplot(alpha = 0.8) +
  scale_x_discrete(labels = c('TRUE' = 'Blackout',
                              'FALSE' = 'No blackout')) +
  scale_fill_manual(values = c('TRUE' = 'firebrick',
                               'FALSE' = 'midnightblue'),
                    labels = c('Blackout', 'No blackout')) +
  labs(title = 'Median household income for census tracts',
       subtitle = 'Blackout status in Houston census tracts',
       x = 'Blackout status',
       y = 'Median household income in $',
       fill = 'Status') +
  scale_y_continuous(labels = scales::dollar_format()) +
  theme_minimal() +
  theme(plot.title = element_text(face = 'bold', size = 16),
        plot.subtitle = element_text(size = 14),
        legend.position = 'bottom')

Figure 1: This plot compares the distributions of median household income between census tracts that experienced blackouts and those that did not during the severe winter storms of February 2021 in Houston, Texas.

Reflection on blackout impacts and study limitations

This analysis used satellite data to compare median household incomes between areas affected and unaffected by the blackout. Interestingly, higher-income census tracts were more affected, as their median income was slightly higher than that of non-blackout areas. This result may reflect several limitations such as wealthier neighborhoods often have more outdoor lighting, census tract–level data may overlook neighborhood differences, and blackout areas could be misclassified due to highway lighting or cloud cover. The analysis may also have assumed all structures were residential, while vegetation near power lines and the 200-meter buffer used could further limit accuracy of true lighting.